rm(list = ls())
library(dplyr)
library(texreg)
library(survival)


#remove everything
rm(list = ls())

# Data ####
d <- read.csv("data.csv")
d <- d %>% dplyr::select(SideA, SideB, DyadId, ConflictId, GWNoA, Year, 
                         Negotiation, Negotiation_lag, Negotiation_inaugural,
                         demref, 
                         polyarchy, duration_dy_ln, imp.gdppcln, imp.popln,
                         Intensity, 
                         Region, RebStrBin.imp, Ethnic.incomp, Incompatibility, 
                         medprev2y,
                         ExChange, oda.pc.ln, UNPKO.total.ln, BdBest_ln, tsln, tsln2, tsln3,
                         demref_poly_xp2_ip2_lag,
                         demref_lexical_autlib_dec5_lag, demref_lexical_demtrans_dec5_lag)

# Factorize
d$Intensity <- as.factor(d$Intensity)
d$Region <- as.factor(d$Region)
d$Incompatibility <- as.factor(d$Incompatibility)

# Figure 3: Mean first differences ####



# Robustness: Replication with Lexical ####

# Model outputs are reported in Appendix as well
# Figure 7 is done using Lexical dataset
# Alternative measurement: authoritarian liberalization and democratic transition
# Main text explains the logic behind this alternative measurement
# We replicate the results with this alternative measurement

# authoritarian liberalization
d$autlib   <- d$demref_lexical_autlib_dec5_lag
# democratic transition
d$demtrans <- d$demref_lexical_demtrans_dec5_lag

# it is possible to have both due to decay function. take only the most recent
d$demtrans[d$autlib > d$demtrans] <- 0
d$autlib[d$autlib < d$demtrans] <- 0

# Formulas:
#model 1 
fx1 <- as.formula (Negotiation ~  autlib + demtrans + 
                      duration_dy_ln  + imp.gdppcln + imp.popln + 
                      Intensity + Region + tsln + tsln2 + tsln3)
#model 2
fx2 <- as.formula (Negotiation ~   autlib + demtrans + 
                      duration_dy_ln  +  imp.gdppcln + imp.popln + 
                      Intensity + Region + 
                      RebStrBin.imp + Ethnic.incomp  + medprev2y +
                      tsln + tsln2 + tsln3  )

#model 3
fx3 <- as.formula (Negotiation ~ duration_dy_ln  + imp.gdppcln + imp.popln + Intensity +
                      Region + RebStrBin.imp + Ethnic.incomp  + oda.pc.ln + medprev2y  +
                      tsln + tsln2 + tsln3 + BdBest_ln + UNPKO.total.ln + autlib + demtrans)

#model 4: fixed effect
fx4 <- as.formula(Negotiation ~ autlib + demtrans +
                     duration_dy_ln  + imp.gdppcln + imp.popln + Intensity +
                     medprev2y  + tsln + tsln2 + tsln3 + strata(DyadId))

# Models 7-8  Transition
fx7 <- as.formula (Negotiation ~  autlib + demtrans +
                      duration_dy_ln  +  imp.gdppcln + imp.popln + Intensity +
                      Region + RebStrBin.imp + Ethnic.incomp + medprev2y )

# Model 9: inaugural same with 7-8
fx9 <- as.formula (Negotiation_inaugural ~  autlib + demtrans +
                          duration_dy_ln  +  imp.gdppcln + imp.popln + Intensity +
                          Region + RebStrBin.imp + Ethnic.incomp + medprev2y )


# Models
# First 4
mx1 <- glm(formula = fx1, data = d, family = binomial(link = "logit"))
mx2 <- glm(formula = fx2, data = d, family = binomial(link = "logit"))
mx3 <- glm(formula = fx3, data = d, family = binomial(link = "logit"))
mx4 <- survival::clogit(fx4, data = d)

screenreg(list(mx1, mx2, mx3, mx4))

# markov transition

mx7 <- glm(formula = fx7, data = d[d$Negotiation_lag == 0,], family = binomial(link = "logit"))
mx8 <- glm(formula = fx7, data = d[d$Negotiation_lag == 1,], family = binomial(link = "logit"))
screenreg(list(mx7, mx8))

## model 9: inaugural
mx9 <- glm(formula = fx9, data = d[d$Negotiation_inaugural >= 0,], 
           family = binomial(link = "logit"))

screenreg(list(mx7, mx8, mx9))

## decay function for the graph #####
decay <- function(x, halflife){ 
  lambda <- log(2) / halflife
  n  <- length(x) ## length of the vector
  tg <- 1:n ## tg of the t global. from 1 to n, length of the vector. this is the row number for the final output
  m <- matrix(nrow = length(x), ncol = length(x)) ## a square empty matrix. 
  ## loop over rows (each observation in vector x) 
  for(i in tg){
    te <- 1:(n - i + 1) # time for observation, will get into the decay function. iow: t of the decay 
    y  <- x[i] * exp(-lambda * te ) # this is the decay function N_t(t) = N_0  * e^(-tau = t) where N_0 is the initial value
    y  <- c(x[i], y)
    y  <- y[1:(length(y)-1)]
    y  <- round(y, 4)
    m[i:n,i] <- y
    
  }
  output <- apply(m, 1, max, na.rm = T) 
  output[output == -Inf] <- NA
  return(output)
  
}

## Figure 7 ####
n_year <- 5
x2 <- decay(c(1,rep(0,n_year)), 5)

### model matrix
m <- mx1
x <- model.matrix(m)

# baseline: no reform:
# take the column nos 
i_demtrans <- match("demtrans", colnames(x))
i_autlib   <- match("autlib", colnames(x))

# dem trans should be 0 and aut lob should be 0 
x[, i_demtrans] <- 0
x[, i_autlib] <- 0

# array for democratic transition
a  <- array(x, c(dim(x) , length(x2)))
dim(a)


# array for authoritarian liberlization
a2 <- array(x, c(dim(x) , length(x2)))

# fill inside
for(i in  1:dim(a)[3]){
  a[,i_demtrans,i] <- x2[i]
  a2[,i_autlib,i] <- x2[i]
}

## coefdraws
nsim <- 10 # make it 10k

coef.draws <- MASS::mvrnorm(n = nsim, mu = m$coefficients, Sigma = vcov(m))

#dim(coef.draws)

## empty ystar
# no ref
da0_ystar <- matrix(NA, nrow(x), nrow(coef.draws ) )
#dim(da0_ystar)
# dem transition
da1_ystar <- array(NA, dim = c(nrow(x), nrow(coef.draws), dim(a)[3]))
#dim(da1_ystar)
#dim(da0_ystar)
#dim(coef.draws)
# autlib
da2_ystar <- array(NA, dim = c(nrow(x), nrow(coef.draws), dim(a)[3]))

## fill inside of ystars
#no ref
da0_ystar <- x %*% t(coef.draws)
da0_yhat <- exp(da0_ystar) / (exp(da0_ystar) + 1)
rm(da0_ystar)
range(da0_yhat)

## ref
for(i in 1:dim(a)[3]){
  da1_ystar[,,i] <- a[,,i] %*% t(coef.draws)
}
da1_yhat <- exp(da1_ystar) / (exp(da1_ystar) + 1)
rm(da1_ystar)
range(da1_yhat)

## autlib 
for(i in 1:dim(a2)[3]){
  da2_ystar[,,i] <- a2[,,i] %*% t(coef.draws)
}
da2_yhat <- exp(da2_ystar) / (exp(da2_ystar) + 1)
rm(da2_ystar)
range(da2_yhat)

### differences : empty arrays
# differences between no ref and categories of interest(i.e., dem transition and authoritarian liberalization)
yhat_fark <- array(NA, dim(da1_yhat))
yhat_fark2 <- array(NA, dim(da1_yhat))

dim(yhat_fark)
dim(da0_yhat)

## loop over to calculate the differences
# for democratic transition
for(i in 1:dim(a)[3]){
  yhat_fark[,,i] <- da1_yhat[,,i] - da0_yhat
}

## for authoritarian liberalization
for(i in 1:dim(a2)[3]){
  yhat_fark2[,,i] <- da2_yhat[,,i] - da0_yhat
}

## average (mean) mean marginal effect
xx  <- apply(yhat_fark,c(2,3), mean)
xx2 <- apply(yhat_fark2,c(2,3), mean)
dim(xx)
dim(xx2)

## take the quantiles
pe <- apply(xx,2,quantile, 0.50)
lb <- apply(xx,2,quantile, 0.95)
ub <- apply(xx,2,quantile, 0.05)

pe2 <- apply(xx2,2,quantile, 0.50)
lb2 <- apply(xx2,2,quantile, 0.95)
ub2 <- apply(xx2,2,quantile, 0.05)

## lets plot
png(filename = "lexical.png", width = 20, height = 10, res = 300, units = "cm") 
par(mar = c(3,3,2,2), mgp = c(1.5,0.5,0), oma = c(0,0,0,0))
plot(pe, ylim = range(c(0, lb, ub, lb2, ub2, 0.2) ), t = "l", lwd = 3, axes  = F,
     xlab = "Years after Democratic Reform",
     ylab = expression(paste(Delta, "Pr(Negotiation)", sep = "")), xlim = c(1,5),
     main = "Democratic Reform and Negotiations", )
axis(1)
axis(2, at = seq(0,0.20,0.05))
polygon(x = c(1:6,6:1) , y = c(lb[1:6],rev(ub[1:6])), col = "grey", border = "NA") 
lines(pe, t = "l", lwd = 3 , lty = 1)
lines(pe2, t = "l", lwd = 3 , lty = 2, col = "black")
lines(lb2, t = "l", lwd = 1 , lty = 3, col = "black")
lines(ub2, t = "l", lwd = 1 , lty = 3, col = "black")
text(y = 0.19, x = 4.25, "Democratic Transition", adj = 0, cex = 0.75)
polygon( x= c(4,4.2,4.2,4), y = c(0.18, 0.18, 0.2, 0.2), col = "grey", border = "NA") 
segments(y0 = 0.19, y1 = 0.19, x0 = 4, x1 = 4.2, lwd =3)
segments(y0 = 0.17, y1 = 0.17, x0 = 4, x1 = 4.2, lwd =1, lty = 3)
segments(y0 = 0.16, y1 = 0.16, x0 = 4, x1 = 4.2, lwd =2, lty = 2)
segments(y0 = 0.15, y1 = 0.15, x0 = 4, x1 = 4.2, lwd =1, lty= 3)
text(y = 0.16, x = 4.25, "Autocratic Liberalization", adj = 0, cex = 0.75)
dev.off()